##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(RColorBrewer)
library(ggpubr)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--cnv_file"), type = "character") ,
    make_option(c("--igc_maf_file"), type = "character") ,
    make_option(c("--dgc_maf_file"), type = "character") ,
    make_option(c("--njmu_info"), type = "character") ,
    make_option(c("--public_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    cnv_file <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/revise/TP53/TP53.cnv.bed"
    njmu_info <- "~/20220915_gastric_multiple/dna_combinePublic/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.rmMIX.tsv"
    public_info <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/combine/MutationInfo.combine.addMolecularSubType.rmMIX.tsv"
    igc_maf_file <- "~/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.IGC.maf"
    dgc_maf_file <- "~/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.DGC.maf"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/revise/TP53"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

njmu_info <- opt$njmu_info
public_info <- opt$public_info
igc_maf_file <- opt$igc_maf_file
dgc_maf_file <- opt$dgc_maf_file
out_path <- opt$out_path
cnv_file <- opt$cnv_file

###########################################################################################

dir.create(out_path , recursive = T)
col <- c( "#006699","#DDA520"  )

###########################################################################################

dat_info_njmu <- data.frame(fread(njmu_info))
dat_info_public <- data.frame(fread(public_info))
dat_cnv <- data.frame(fread(cnv_file))
dat_igc <- data.frame(fread(igc_maf_file))
dat_dgc <- data.frame(fread(dgc_maf_file))

###########################################################################################

Variant_Types <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")
dat_igc <- subset( dat_igc , Variant_Classification %in% Variant_Types & Hugo_Symbol == "TP53" )
dat_dgc <- subset( dat_dgc , Variant_Classification %in% Variant_Types & Hugo_Symbol == "TP53" )

mut_sample <- c(dat_igc$Tumor_Sample_Barcode , dat_dgc$Tumor_Sample_Barcode)
msi_sample <- subset( dat_info_public , MS_Type == "MSI" )$Tumor

###########################################################################################

dat_info_njmu$id <- paste0(dat_info_njmu$Tumor , "_" , dat_info_njmu$Normal)
dat_cnv_njmu <- merge( dat_cnv , dat_info_njmu[,c("id" , "Patient" , "Class")] , by.x = "Sample" , by.y = "id" )
dat_info_public$id <- dat_info_public$Tumor
dat_info_public$Patient <- dat_info_public$Tumor
dat_cnv_tcga <- merge( dat_cnv , dat_info_public[,c("id" , "Patient" , "Class")] , by.x = "Sample" , by.y = "id" )
dat_cnv <- rbind(dat_cnv_njmu , dat_cnv_tcga)

###########################################################################################
result_plot <- c()
for(ClassN in unique(dat_cnv$Class)){
    print(ClassN)
    
    tmp <- subset( dat_cnv , Class == ClassN )
   
    ## 判断拷贝数改变
    tmp$CNV_type <- "Nochange"
    tmp$CNV_type <- ifelse( tmp$Corrected_MinorCN == 0 & tmp$Corrected_Copy_Number==2 , "LOH" , tmp$CNV_type )
    tmp$CNV_type <- ifelse( tmp$Corrected_Copy_Number > 2 , "Gain" , tmp$CNV_type )
    tmp$CNV_type <- ifelse( tmp$Corrected_Copy_Number <= 1 , "Loss" , tmp$CNV_type )

    res_tmp <- c()
    for( sample in unique(tmp$Patient) ) {
        tmp_use <- subset( tmp , Patient == sample )
        tmp_use <- unique(tmp_use[,c("Patient" , "CNV_type")])
        ## 三个样本存在多个患者拷贝数改变不一致，选择有改变的样本
        if( nrow(tmp_use) > 1 ){
            tmp_use <- tmp_use[1,]
        }
        res_tmp <- rbind(res_tmp , tmp_use)
    }

    res_tmp$CNV_change <- ifelse( res_tmp$CNV_type == "Nochange" , "Nochange" , "DoubleHit" )
    res_tmp$Class <- ClassN

    result_plot <- rbind( result_plot , res_tmp )
}

result_plot$mut <- ifelse( result_plot$Patient %in% mut_sample , "Mut" , "Wild" )
result_plot <- subset( result_plot, !Patient %in% msi_sample )

out_name <- paste0(out_path , "/Bialleic.TP53.csv")
write.csv( result_plot , out_name , row.names = F , quote = F )
result_plot$bialleic <- ifelse( result_plot$CNV_type %in% c("LOH" , "Loss") & result_plot$mut == "Mut" , "Bialleic\n(Mutation with LOH or Loss)" , "Other" )

###########################################################################################
## 基因分布堆叠图
if(1!=1){
    res_tmp2 <- data.frame(table(result_plot$Class , result_plot$CNV_type , result_plot$mut))
    colnames(res_tmp2) <- c("Class" , "CNV_type" , "Mut_type" , "count")
    res_tmp2 <- res_tmp2 %>%
    group_by(Class) %>%
    mutate(count_all=sum(count),
       ratio=count/count_all)
    res_tmp2 <- subset( res_tmp2 , Class != "IM" )
}

res_tmp2 <- data.frame(table(result_plot$Class , result_plot$bialleic))
colnames(res_tmp2) <- c("Class" , "bialleic" , "count")
res_tmp2 <- res_tmp2 %>%
group_by(Class) %>%
mutate(count_all=sum(count),
   ratio=count/count_all)
res_tmp2 <- subset( res_tmp2 , Class != "IM" )

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("P == ",down," %*% 10","^",up)
    return(text)
}

p <- fisher.test(matrix(as.numeric(res_tmp2[order(res_tmp2$Class),]$count) , ncol =2))$p.value
if( p < 0.001 ){
    p_text <- trans(p)
}else{
    p_text <- paste0( "P == " , round(as.numeric(p) , 3) ) 
}

res_tmp2$p <- ""
res_tmp2$p[1] <- p
res_tmp2$p_text <- ""
res_tmp2$p_text[1] <- p_text

result <- data.frame(res_tmp2)
result$ratio[is.na(result$ratio)] <- 0
result$value_text <- round(result$ratio , 2) * 100

#result$CNV_type <- factor( result$CNV_type , levels = c("LOH" , "Loss" , "Gain" , "Nochange") )
#result$Mut_type <- factor( result$Mut_type , levels = c("Mut" , "Wild") )
result$Class <- factor( result$Class , levels = c("IGC" , "DGC") )

col <- c(
    rgb(190,97,99,alpha=255,maxColorValue=255),
    rgb(91,139,101,alpha=255,maxColorValue=255),
    rgb(76,123,161,alpha=255,maxColorValue=255),
    "grey"
  )

names(col) <- c("Gain" , "LOH" , "Loss" , "Nochange")
col <- c(
    rgb(190,97,99,alpha=255,maxColorValue=255),
    "grey"
  )

names(col) <- c("Bialleic\n(Mutation with LOH or Loss)" , "Other")
result$value_text <- ifelse( result$ratio == 0 , "" , result$value_text )

plot <- ggplot( data = result , aes( x = Class , y = ratio , fill = bialleic ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion (%)")+
    theme(panel.grid = element_blank())+
    scale_fill_manual(values=col) +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 8,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold') ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/Bialleic.TP53.pdf")
ggsave( out_name , plot , width = 4.2 , height = 4 )
out_name <- paste0(out_path , "/Bialleic.TP53.p.csv")
write.csv( result , out_name , row.names = F , quote = F )